data <- readRDS("~/Downloads/week_11/data/historical_web_data_26112015.rds")
plot_ly(data, x = data$Latitude, y = data$Longitude, z = data$Depth, size = data$Magnitude) %>% layout(scene = list(xaxis = list(title = 'Latitude'), yaxis = list(title = 'Longitude'), zaxis = list(title = 'Depth')))
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter3d
## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
bigEQ = read_delim('~/Downloads/week_11/data/disaster.txt', '\t', escape_double = FALSE, trim_ws = TRUE)
sunami = bigEQ %>% filter(FLAG_TSUNAMI == 'Tsu')
sunami = sunami[!is.na(sunami$EQ_PRIMARY),]
sunami = sunami %>% filter(EQ_PRIMARY >= 8)
worldMap = map_data('world')
ggplot() + geom_polygon(data = worldMap, aes(x = long, y = lat, group = group),fill = 'blue') + geom_point(data = sunami, aes(x = LONGITUDE, y = LATITUDE), color = 'orange', size = 0.3) + transition_time(EQ_PRIMARY) + labs(title = 'mag {frame_time}') + ease_aes('linear')
iranEQ <- readRDS("~/Downloads/week_11/data/iran_earthquake.rds")
ggmap(read_rds("~/Downloads/week_11/data/Tehrn_map_6.rds")) + stat_density_2d(geom = 'polygon', data = iranEQ, aes(x = Long, y = Lat, fill = ..level.., alpha = ..level..))
## Warning: Removed 243 rows containing non-finite values (stat_density2d).
زلزله بزرگ را زلزله با قدرت بیشتر از ۶.۷ ریشتر در نظر گرفتیم.
irBig = bigEQ %>% filter(EQ_PRIMARY > 6 & COUNTRY == 'IRAN')
distanceYear = data.frame(irBig[-1,]$YEAR - irBig[-nrow(irBig),]$YEAR)
colnames(distanceYear) = c('year')
notHappened = distanceYear %>% filter(year > 5)
cat(1 - nrow(notHappened)/nrow(distanceYear))
## 0.6764706
meannumeq <- bigEQ %>% group_by(COUNTRY) %>% summarise(numberOfEQ = n(), meanDeath = sum(DEATHS,na.rm = T)/numberOfEQ, alldeath = sum(DEATHS,na.rm = T))
mapDevice('x11')
heatMap = joinCountryData2Map(meannumeq, joinCode="NAME", nameJoinColumn="COUNTRY")
## 138 codes from your data successfully matched countries in the map
## 16 codes from your data failed to match with a country code in the map
## 105 codes from the map weren't represented in your data
mapCountryData(heatMap, nameColumnToPlot="alldeath")
## You asked for 7 quantiles, only 5 could be created in quantiles classification
fittedLine <- lm(bigEQ, formula = DEATHS ~ LONGITUDE + LATITUDE + FOCAL_DEPTH + EQ_PRIMARY )
summary(fittedLine)
##
## Call:
## lm(formula = DEATHS ~ LONGITUDE + LATITUDE + FOCAL_DEPTH + EQ_PRIMARY,
## data = bigEQ)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12556 -3519 -1872 35 311710
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -15451.303 3081.358 -5.014 6.13e-07 ***
## LONGITUDE -1.232 5.813 -0.212 0.83226
## LATITUDE 64.237 21.852 2.940 0.00335 **
## FOCAL_DEPTH -21.929 11.354 -1.931 0.05367 .
## EQ_PRIMARY 2678.740 464.175 5.771 1.01e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15440 on 1178 degrees of freedom
## (4843 observations deleted due to missingness)
## Multiple R-squared: 0.031, Adjusted R-squared: 0.0277
## F-statistic: 9.42 on 4 and 1178 DF, p-value: 1.7e-07
worldEQ = read.csv('~/Downloads/week_11/data/worldwide.csv')
worldEQ$date = as_date(worldEQ$time)
worldEQ$ID = interaction(day(worldEQ$date), year(worldEQ$date), month(worldEQ$date), sep='')
worldEQ$country = sub(".*, ", "", worldEQ$place)
left = worldEQ %>% group_by(country, ID) %>% filter(mag == max(mag))
right = worldEQ %>% group_by(country, ID) %>% mutate(ind = rank(desc(mag))) %>% arrange(ind) %>% filter(ind == 2)
inter = inner_join(left, right, by = c('country', 'ID'))
train = inter[1:as.integer(0.9 * nrow(inter)),]
test = inter[-(1:as.integer(0.9 * nrow(inter))),]
model = lm(data = train, mag.x ~ mag.y)
summary(model)
##
## Call:
## lm(formula = mag.x ~ mag.y, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.39643 -0.28619 -0.09373 0.11327 2.61112
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.369494 0.022151 16.68 <2e-16 ***
## mag.y 1.005386 0.005515 182.31 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3763 on 7792 degrees of freedom
## Multiple R-squared: 0.8101, Adjusted R-squared: 0.8101
## F-statistic: 3.324e+04 on 1 and 7792 DF, p-value: < 2.2e-16
summary(stats::predict(model, test) - test$mag.x)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.2057 -0.1090 0.1338 0.0153 0.2921 0.3943
According to correlation coefficient and scatter plot, there is not any relation between depth and magnitude of earthquak
cor.test(worldEQ$depth,worldEQ$mag)
##
## Pearson's product-moment correlation
##
## data: worldEQ$depth and worldEQ$mag
## t = 50.03, df = 60629, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1914584 0.2067469
## sample estimates:
## cor
## 0.1991147
worldEQ = read.csv('~/Downloads/week_11/data/worldwide.csv')
worldEQ$country = sub(".*, ", "", worldEQ$place)
worldEQ$time = as.Date(worldEQ$time)
worldEQ$year = year(worldEQ$time)
worldEQ = worldEQ %>% group_by(place, year) %>% summarise(number = n()) %>% ungroup() %>% group_by(place) %>% summarise(n = mean(number))
worldEQ %>% arrange(desc(n)) %>% top_n(20) -> top20
## Selecting by n
ggplot(top20) + geom_bar(aes(x = place, y = n), stat = "identity")
bigEQ %>% arrange(desc(EQ_PRIMARY)) %>% top_n(10) ->f
## Selecting by TOTAL_HOUSES_DAMAGED_DESCRIPTION
ggplot(f) + geom_bar(aes(x = COUNTRY, y = EQ_PRIMARY),stat = "identity")
## Warning: Removed 3 rows containing missing values (position_stack).
bigEQ %>% filter(EQ_PRIMARY > 7) %>% group_by(COUNTRY) %>% summarise(n = n()) %>% top_n(10) -> d
## Selecting by n
ggplot(d) + geom_bar(aes(x = COUNTRY, y = n),stat = "identity")
ts = sunami %>% filter(EQ_PRIMARY > 8)
ggplot() + geom_polygon(data = worldMap, aes(x = long, y = lat, group = group),fill = 'blue') + geom_point(data = ts, aes(x = LONGITUDE, y = LATITUDE), color = 'orange', size = 0.3)